%% 多名运动员蹦床疲劳损伤分析与优化
% 基于题目要求构建模型，分析多人蹦床运动时的受力情况和疲劳损伤，
% 并提出优化策略以减少蹦床疲劳损伤

clear all; close all; clc;

%% 定义参数
% 运动员参数
athletes = struct();
athletes(1).gender = 'male';   athletes(1).height = 1.75; athletes(1).leg = 0.88; athletes(1).weight = 68; athletes(1).pos = [2.0, 1.0]; athletes(1).COM_height = 0.93; athletes(1).k = 4951.05;
athletes(2).gender = 'female'; athletes(2).height = 1.63; athletes(2).leg = 0.81; athletes(2).weight = 55; athletes(2).pos = [3.8, 0.5]; athletes(2).COM_height = 0.84; athletes(2).k = 4932.36;
athletes(3).gender = 'male';   athletes(3).height = 1.80; athletes(3).leg = 0.90; athletes(3).weight = 72; athletes(3).pos = [1.2, 1.5]; athletes(3).COM_height = 0.95; athletes(3).k = 4970.56;
athletes(4).gender = 'female'; athletes(4).height = 1.68; athletes(4).leg = 0.84; athletes(4).weight = 60; athletes(4).pos = [0.8, 1.8]; athletes(4).COM_height = 0.86; athletes(4).k = 4973.62;
athletes(5).gender = 'male';   athletes(5).height = 1.70; athletes(5).leg = 0.85; athletes(5).weight = 65; athletes(5).pos = [3.0, 1.2]; athletes(5).COM_height = 0.88; athletes(5).k = 4935.29;

% 蹦床参数
trampoline.length = 4.028;    % 蹦床网长度(m)
trampoline.width = 2.014;     % 蹦床网宽度(m)
trampoline.springs = 112;     % 弹簧数量
trampoline.g = 9.81;          % 重力加速度(m/s^2)
trampoline.damping = 0.2;     % 阻尼系数(假设值)
trampoline.max_deflection = 1.0; % 最大变形量(m)
trampoline.restitution = 0.8; % 恢复系数(假设值)

% 疲劳模型参数 (基于修正的Miner准则)
trampoline.fatigue.stress_limit = 500;      % 疲劳极限应力(N/m^2)
trampoline.fatigue.C = 1e6;                 % 疲劳寿命系数
trampoline.fatigue.m = 3;                   % 疲劳指数(越大疲劳越快)
trampoline.fatigue.initial_life = 1;        % 初始疲劳寿命(归一化)

% 仿真参数
sim.dt = 0.001;                  % 时间步长(s)
sim.total_time = 30;             % 总仿真时间(s)
sim.jump_height = 3;             % 标准起跳高度(m)
sim.min_jump_interval = 0.5;     % 最小起跳间隔(s)

%% 1. 建立动力学方程
% 定义蹦床变形计算函数
function [deflection, force] = calculate_deflection(athletes, positions, trampoline, dt)
    n = length(athletes);
    deflection = zeros(n, 1);
    force = zeros(n, 1);
    
    for i = 1:n
        if positions(i, 3) <= 0
            % 计算接触点的垂直形变
            deflection(i) = -positions(i, 3);
            
            % 计算弹力和阻尼力
            spring_force = athletes(i).k * deflection(i);
            damping_force = trampoline.damping * (positions(i, 6) < 0) * abs(positions(i, 6));
            
            % 总力
            force(i) = spring_force + damping_force;
        end
    end
end

% 定义状态更新函数
function new_positions = update_positions(athletes, positions, forces, trampoline, dt)
    n = length(athletes);
    new_positions = positions;
    
    for i = 1:n
        % 更新状态 [x, y, z, vx, vy, vz]
        % 计算合力
        Fx = 0;  % 水平方向无力
        Fy = 0;  % 水平方向无力
        
        if positions(i, 3) <= 0
            % 在蹦床上
            Fz = forces(i) - athletes(i).weight * trampoline.g;
        else
            % 空中
            Fz = -athletes(i).weight * trampoline.g;
        end
        
        % 更新加速度
        ax = Fx / athletes(i).weight;
        ay = Fy / athletes(i).weight;
        az = Fz / athletes(i).weight;
        
        % 更新速度
        new_positions(i, 4) = positions(i, 4) + ax * dt;
        new_positions(i, 5) = positions(i, 5) + ay * dt;
        new_positions(i, 6) = positions(i, 6) + az * dt;
        
        % 更新位置
        new_positions(i, 1) = positions(i, 1) + new_positions(i, 4) * dt;
        new_positions(i, 2) = positions(i, 2) + new_positions(i, 5) * dt;
        new_positions(i, 3) = positions(i, 3) + new_positions(i, 6) * dt;
    end
end

%% 2. 分析不同条件下的蹦床受力和疲劳损伤

% 定义疲劳损伤计算函数
function damage = calculate_fatigue_damage(max_forces, trampoline)
    % 基于修正的Miner准则计算疲劳损伤
    % 每次循环的损伤 = 1/N，其中N是在该应力下的疲劳寿命
    
    damage = 0;
    for i = 1:length(max_forces)
        if max_forces(i) > 0
            % 计算等效应力
            stress = max_forces(i) / (trampoline.length * trampoline.width);
            
            % 计算该应力水平下的疲劳寿命
            if stress > trampoline.fatigue.stress_limit
                N = trampoline.fatigue.C / (stress^trampoline.fatigue.m);
                
                % 累积损伤
                damage = damage + 1/N;
            end
        end
    end
end

% 定义仿真函数
function [max_forces, total_damage] = simulate_trampoline(athletes, jump_times, trampoline, sim)
    n = length(athletes);
    time_steps = round(sim.total_time / sim.dt);
    
    % 初始位置和速度 [x, y, z, vx, vy, vz]
    positions = zeros(n, 6);
    for i = 1:n
        positions(i, 1:2) = athletes(i).pos;  % x,y坐标
        positions(i, 3) = athletes(i).COM_height;  % 初始高度
    end
    
    % 记录每个时间步的最大力
    forces_history = zeros(time_steps, 1);
    
    % 模拟每个时间步
    for t = 1:time_steps
        current_time = t * sim.dt;
        
        % 检查是否有运动员在此时间点起跳
        for i = 1:n
            if abs(current_time - jump_times(i)) < sim.dt && positions(i, 3) <= athletes(i).COM_height
                % 设置起跳初速度
                jump_velocity = sqrt(2 * trampoline.g * sim.jump_height);
                positions(i, 6) = jump_velocity;
            end
        end
        
        % 计算蹦床变形和力
        [deflections, forces] = calculate_deflection(athletes, positions, trampoline, sim.dt);
        
        % 更新位置和速度
        positions = update_positions(athletes, positions, forces, trampoline, sim.dt);
        
        % 记录当前时间步的总力
        forces_history(t) = sum(forces);
    end
    
    % 计算最大力
    max_forces = forces_history;
    
    % 计算总疲劳损伤
    total_damage = calculate_fatigue_damage(max_forces, trampoline);
end

%% 3. 提出优化策略
% 基于遗传算法优化起跳时序

% 定义适应度函数
function fitness = evaluate_strategy(jump_times, athletes, trampoline, sim)
    [~, total_damage] = simulate_trampoline(athletes, jump_times, trampoline, sim);
    fitness = -total_damage;  % 最小化疲劳损伤
end

% 定义遗传算法优化函数
function [best_jump_times, best_fitness] = optimize_jump_strategy(athletes, trampoline, sim)
    n = length(athletes);
    
    % 遗传算法参数
    pop_size = 50;
    max_gen = 100;
    mutation_rate = 0.1;
    crossover_rate = 0.8;
    
    % 初始化种群
    population = zeros(pop_size, n);
    for i = 1:pop_size
        % 随机产生起跳时间序列(所有人必须在出现落地者前起跳)
        % 假设一个人从起跳到落地的时间约为1秒
        max_jump_time = 2;  % 最晚的起跳时间
        
        % 随机排列起跳顺序
        jump_order = randperm(n);
        
        % 设定起跳时间
        for j = 1:n
            population(i, jump_order(j)) = (j-1) * sim.min_jump_interval + rand() * sim.min_jump_interval;
        end
    end
    
    % 初始化适应度
    fitness = zeros(pop_size, 1);
    
    % 开始优化
    for gen = 1:max_gen
        % 评估适应度
        for i = 1:pop_size
            fitness(i) = evaluate_strategy(population(i, :), athletes, trampoline, sim);
        end
        
        % 记录最优解
        [best_fitness, best_idx] = max(fitness);
        best_jump_times = population(best_idx, :);
        
        % 创建新一代
        new_population = zeros(pop_size, n);
        
        % 精英保留
        elite_count = 2;
        [~, elite_indices] = sort(fitness, 'descend');
        new_population(1:elite_count, :) = population(elite_indices(1:elite_count), :);
        
        % 生成新个体
        for i = elite_count+1:pop_size
            % 选择父母
            parent1_idx = tournament_selection(fitness);
            parent2_idx = tournament_selection(fitness);
            
            % 交叉
            if rand() < crossover_rate
                [child1, ~] = crossover(population(parent1_idx, :), population(parent2_idx, :));
                new_population(i, :) = child1;
            else
                new_population(i, :) = population(parent1_idx, :);
            end
            
            % 变异
            new_population(i, :) = mutate(new_population(i, :), mutation_rate, sim.min_jump_interval, max_jump_time);
        end
        
        % 更新种群
        population = new_population;
        
        % 输出进度
        if mod(gen, 10) == 0
            fprintf('Generation %d: Best fitness = %f\n', gen, best_fitness);
        end
    end
end

% 定义锦标赛选择函数
function selected_idx = tournament_selection(fitness)
    k = 3;  % 锦标赛大小
    indices = randi(length(fitness), 1, k);
    [~, best_of_k] = max(fitness(indices));
    selected_idx = indices(best_of_k);
end

% 定义交叉函数
function [child1, child2] = crossover(parent1, parent2)
    n = length(parent1);
    crossover_point = randi(n-1);
    
    child1 = [parent1(1:crossover_point), parent2(crossover_point+1:end)];
    child2 = [parent2(1:crossover_point), parent1(crossover_point+1:end)];
end

% 定义变异函数
function child = mutate(child, mutation_rate, min_interval, max_time)
    n = length(child);
    
    for i = 1:n
        if rand() < mutation_rate
            % 在允许范围内变异
            child(i) = rand() * max_time;
        end
    end
    
    % 确保满足起跳时间最小间隔约束
    child = sort(child);
    for i = 2:n
        if child(i) - child(i-1) < min_interval
            child(i) = child(i-1) + min_interval;
        end
    end
end

%% 4. 主函数：评估策略效果

% 定义默认起跳时序(均匀分布)
default_jump_times = zeros(1, length(athletes));
for i = 1:length(athletes)
    default_jump_times(i) = (i-1) * 2;  % 每2秒一个人起跳
end

% 模拟默认策略
[default_forces, default_damage] = simulate_trampoline(athletes, default_jump_times, trampoline, sim);
fprintf('默认策略疲劳损伤: %f\n', default_damage);

% 优化起跳时序
[optimized_jump_times, best_fitness] = optimize_jump_strategy(athletes, trampoline, sim);
fprintf('优化策略疲劳损伤: %f\n', -best_fitness);

% 模拟优化策略
[optimized_forces, optimized_damage] = simulate_trampoline(athletes, optimized_jump_times, trampoline, sim);

% 计算寿命提升效果
life_improvement = default_damage / optimized_damage;
fprintf('寿命提升比例: %.2f倍\n', life_improvement);

% 可视化结果
figure;
subplot(2,1,1);
t = (1:length(default_forces)) * sim.dt;
plot(t, default_forces);
title('默认策略下蹦床受力');
xlabel('时间(s)');
ylabel('力(N)');

subplot(2,1,2);
t = (1:length(optimized_forces)) * sim.dt;
plot(t, optimized_forces);
title('优化策略下蹦床受力');
xlabel('时间(s)');
ylabel('力(N)');

figure;
bar([default_damage, optimized_damage]);
set(gca, 'xticklabel', {'默认策略', '优化策略'});
ylabel('疲劳损伤');
title('疲劳损伤比较');

figure;
bar(1:length(athletes), [default_jump_times; optimized_jump_times]);
legend('默认起跳时间', '优化起跳时间');
xlabel('运动员编号');
ylabel('起跳时间(s)');
title('起跳时序比较');

%% 5. 功能函数：蹦床网格绘制与变形可视化
function visualize_trampoline_deformation(athletes, positions, trampoline)
    % 创建网格
    [X, Y] = meshgrid(linspace(0, trampoline.length, 50), linspace(0, trampoline.width, 25));
    Z = zeros(size(X));
    
    % 计算每个网格点的变形
    for i = 1:length(athletes)
        if positions(i, 3) <= 0
            % 计算变形
            for j = 1:size(X, 1)
                for k = 1:size(X, 2)
                    % 计算与运动员位置的距离
                    dist = sqrt((X(j,k) - positions(i,1))^2 + (Y(j,k) - positions(i,2))^2);
                    
                    % 使用高斯函数模拟变形扩散
                    sigma = 0.5;  % 变形范围参数
                    Z(j,k) = Z(j,k) - positions(i,3) * exp(-dist^2/(2*sigma^2));
                end
            end
        end
    end
    
    % 绘制3D表面
    figure;
    surf(X, Y, Z);
    title('蹦床变形可视化');
    xlabel('X(m)');
    ylabel('Y(m)');
    zlabel('Z(m)');
    colorbar;
end

%% 6. 预测蹦床疲劳寿命
function predict_fatigue_life(trampoline, default_damage, optimized_damage)
    % 假设每天训练次数
    jumps_per_day = 1000;
    
    % 计算默认策略下的疲劳寿命(天)
    default_life_days = 1 / (default_damage * jumps_per_day);
    
    % 计算优化策略下的疲劳寿命(天)
    optimized_life_days = 1 / (optimized_damage * jumps_per_day);
    
    % 转换为年
    default_life_years = default_life_days / 365;
    optimized_life_years = optimized_life_days / 365;
    
    fprintf('默认策略下蹦床预计使用寿命: %.2f年\n', default_life_years);
    fprintf('优化策略下蹦床预计使用寿命: %.2f年\n', optimized_life_years);
    fprintf('寿命延长: %.2f年\n', optimized_life_years - default_life_years);
end

% 调用寿命预测函数
predict_fatigue_life(trampoline, default_damage, optimized_damage);